A configuration interaction analysis of exchange in double quantum dots 
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We describe in detail a full configuration interaction (CI) method designed to analyze systems 
of quantum dots. This method is capable of exploring large regions of parameter space, like more 
approximate approaches such as Heitler London and Hund Mulliken, though it is not limited to 
weakly coupled dots. In particular, this method is well-suited to the analysis of solid state quantum- 
dot-based qubits, and we consider the case of a double quantum dot (DQD) singlet-triplet qubit. 
Past analyses have used techniques which are either substantially restricted in the regimes they can 
be used, or device specific and unsuited to exploration of a large regions of parameter space. We 
analyze how the DQD exchange energy, which is central to the operation of qubit rotation gates, 
depends on a generic set of system parameters including magnetic field, DQD detuning, dot size, 
and dot separation. We discuss the implications of these results to the construction of real devices. 
We provide a benchmark of the CI by directly comparing results from the CI method with exact 
results for two electrons in a single parabolic potential (dot). 
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I. INTRODUCTION 

The understanding and analysis of devices consisting 
of multiple coupled quantum dots is becoming increas- 
ingly important. This is particularly evident in the area 
of quantum computing, where there are proposals for 
solid state qubits built from double quantum dots. Such 
qubits must interact with each other (at least pairwise) 
to facilitate the two-qubit gates necessary for universal 
quantum computation. Several standard approaches are 
used to simulate multi-dot systems, ranging from classi- 
cal electronics solvers, which are ideal for large structures 
with many-electron quantum dots, to exact diagonaliza- 
tion techniques such as the full configuration interaction 
method, which arc ideal for small systems of few-electron 
dots. In this work, we are concerned with the latter, and 
specifically with understanding the physics underlying a 
general system of quantum dots. We present here a full 
configuration interaction (CI) method which uses a Gaus- 
sian basis and is well-suited to the many-quantum-dot 
devices used to construct solid state qubits. The method 
differs in small but significant ways from the approaches 
used in previous work, and is applied to gain insight into 
the double quantum dot (DQD) structures relevant to 
quantum computing. It is capable of extending not only 
quantitatively upon prior work limited to more approx- 
imate methods, but also qualitatively by capturing new 
aspects of DQD behavior. In Ref. 1, for instance, the 
(0,2)-occupation regime is studied in detail with respect 
to noise robustness in DQD qubits. The DQD results 
we present here are more than just a slight improvement 
on previous, more approximate calculations; they show 
novel behavior that has not been explored to any great 
extent. Before studying the DQD, we provide a direct 
comparison between the CI and a known exact result 
to help assess the degree of accuracy achieved by the 
method. 

In section II we describe the method in detail, high- 
lighting its differences from similar past approaches. In 
section III, we benchmark the CI using an exactly solv- 



able model: two electrons in a 2D parabolic potential. 
We consider there the convergence of the method with 
basis size, and find that with a relatively small basis quite 
accurate 2-electron energies are obtained. In section IV 
the CI is applied to a singlet-triplet qubit. 2 ' 3 We study 
the exchange energy J as a function of parameters which 
specify the electrostatic potential's shape and the mag- 
netic field. Past studies of the exchange energy in double- 
dot structures have been limited to either more approxi- 
mate methods (e.g. Heitler London (HL) and Hund Mul- 
liken (HM)) or to sophisticated finite-clement numerical 
simulations of specific device structures. Here we are able 
to give a more comprehensive picture of DQD physics 
since the CI is operable in regimes inaccessible to HL 
and HM approaches while still maintaining the compu- 
tational speed necessary for extensive traversal of param- 
eters space. 



II. CALCULATION 

We now describe the configuration interaction (CI) 
method used to solve a ro-electron effective mass Hamil- 
tonian of the form 



(1) 



where k is an effective dielectric constant, and the single 
particle Hamiltonian for the i th particle is given by 
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where r and p are the position and momentum, respec- 
tively, of the i th electron. V is the single-particle poten- 
tial function, and m* is the effective mass (generally a 
tensor). A vector potential A determines the magnetic 
field B = V x A, which we restrict to be constant and 
along the z-direction: B = Bz. 
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Overall, the approach follows the standard CI 
prescription 4 of finding a single-particle basis, building 
many-particle states from this basis, and exactly diago- 
nalizing the Hamiltonian in the resulting many-particle 
basis. In the present approach the single particle ba- 
sis states are linear combinations of uq s-type Gaussian 
functions (in real space) of the form 

g(x, y, z) = Ne- cl *( x - x °) 2 e- a y( y - yo ') 2 , . 

x e -a^(z-z ) 2 e 1 §§-{y x-x y) ^ \ > 

where N is a normalization factor, B is the magnetic 
field, and f = (x ,y , z ) and a = (a x ,a y ,a z ) are the 
position and exponential coefficient of the Gaussian func- 
tion, respectively. The aim in using such functions is 
to reduce the basis size needed to accurately model a 
system while allowing nearly all of the Hamiltonian ma- 
trix elements to be computed analytically, which results 
in a substantial performance advantage over approaches 
which discretize the system using a dense real-space 
mesh. Derivations and final expressions for the relevant 
matrix elements arc given in Appendix A. Once the loca- 
tions and exponential factors of the Gaussians are set (see 
below), an orthogonal single-particle basis is found by di- 
agonalizing the single-particle Hamiltonian (Eq. 2) in the 
non-orthogonal basis of Gaussian functions (a generalized 
eigenvalue problem). All the two-particle Slater determi- 
nant states from the single-particle basis are formed, and 
used as a many-body basis. The full Hamiltonian (Eq. 1) 
is diagonalizcd in this basis, resulting in the system en- 
ergies and eigenstates. The set of Gaussian functions is 
optimized over a subset of fixed number and arrange- 
ment as follows. The number, arrangement, and initial 
locations and widths of the Gaussian functions are in- 
put at the beginning of each run. The set of Gaussian 
functions is then optimized by changing their locations 
and widths (while maintaining the same overall arrange- 
ment) to minimize the many-body ground state energy 
of a given symmetry sector of the Hamiltonian. When 
states with different symmetries are of interest, the min- 
imization is done multiple times (or, equivalently, the 
Hamiltonian matrix could be block-diagonalized and ba- 
sis optimization performed within each block separately) . 
In the case of exchange calculations, where the lowest en- 
ergy singlet and unpolarized triplet states are required, 
the singlet and unpolarized triplet energies are computed 
using separate basis optimizations. Details of this opti- 
mization are given in Appendix B. 



III. SINGLE PARABOLIC DOT 

We first consider an important test case: a single 
parabolic quantum dot in two dimensions. With two 
interacting electrons, the system energies can be found 
exactly (the computation can be done analytically for 
certain dot confinement energies and numerically for any 
set of parameters using a one-dimensional Schrodingcr 



solver). 5 Comparison with this exact solution provides a 
useful benchmark for the CI, and we are able to analyze 
how well a Gaussian basis is able to capture a strongly- 
correlated two-electron state. By varying the basis size, 
information on the convergence of the CI is obtained. As 
a side remark, we note that the CI will always find the 
correct single-electron ground state energy since one of 
the Gaussian basis elements is chosen to be exactly this 
solution. 

Consider a single parabolic dot, given by the potential 

V(r) = \m*ujlr 2 (4) 

where m* is the effective mass and Hujo is the confinement 
energy of the dot. We insert this potential into Eq. 2 and 
then consider the full Hamiltonian given by Eq. 1, where 
n = 2 is the number of electrons (there is no Coulomb 
repulsion term when n = 1). We use GaAs material 
parameters: n = 12.9, and m* = 0.067 m e . 

We solve the full Hamiltonian exactly using the method 
in Ref. 5 to reduce the problem to an ordinary (one- 
dimensional) Schrodingcr equation by switching to center 
of mass coordinates. We then solve this equation using 
the technique prescribed in Ref. 6. The lowest singlet 
(total spin S = 0) and unpolarized triplet (S = 1, S z =0) 
energies obtained by the CI relative to the exact solution 
are shown in Fig. 1 as a function of basis size. As size 
does not uniquely specify a basis, we give the spatial 
arrangement of the basis elements used in Fig. 2 which is 
referenced by the table in in Fig. 1. Overall, we find for 
a range of dot confinement loq and magnetic field B that 
the CI energies converge to within ss 0.5% of their value 
for basis sizes around 10, and to within w 0.05% for basis 
sizes around 20. Thus, for the dot parameters of Fig. 1 
where the energies are of order 10 meV, convergence is 
obtained to within 50 and 5 /ieV for roughly 10 and 20 
basis elements, respectively. This assumes a good basis 
arrangement (cf. Fig. 2), as a poor choice of where to 
place the basis elements (e.g. all in a single line) will 
clearly not produce converged values even for large basis 
sizes. 



IV. DOUBLE PARABOLIC QUANTUM DOT 

We proceed to consider a double quantum dot (DQD) 
system, which cannot be solved or substantially simpli- 
fied analytically. It is directly relevant to quantum com- 
puting as it can be used to implement the fundamental 
unit for quantum computation, a quantum bit (qubit). 
A qubit can be defined as an effective two-level quantum 
mechanical system with a long coherence time that can 
also be controlled and measured. Thus, the state space 
of a qubit, however complex the actual system, can be 
mapped onto that of a single spin-i. Operations on the 
qubit correspond to rotations of the spin-i qubit, and 
are thus referred to as "qubit rotations" . While architec- 
tures have been proposed which implement qubits using 
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FIG. 1: The Cl-computed lowest lying singlet (Eso) and 
triplet (Eto) energies relative to the exact value as a func- 
tion of basis size. The energy, Their difference, the exchange 
energy J = Eto — Eso is also given relative to the exact value. 
The spatial arrangement of the basis elements for each size 
is given in the table to the right, which gives a letter A-D of 
a spatial plot in Fig.2. A trailing Astrix's (*) indicates that 
there are two basis elements lying on top of one another at 
the center of the dot, and the sum (+) of two letters indicates 
that the elements of the referred to spatial plots are combined 
(elements at the same position have different exponential fac- 
tors). Dot parameters hwo = 3meV and B = 0. 
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trons are localized in a double quantum dot. The dots 
are laterally spaced near a material interface, and litho- 
graphically defined gates in a plane parallel to the inter- 
face create and control the electrostatic potential forming 
the quasi-two-dimensional dots. We assume the dot po- 
tentials are parabolic, which is makes the system particu- 
larly suited to analysis by our Gaussian-basis CI method. 
We use this example to both illustrate the capabilities of 
the method described in section II and to report, using 
GaAs material parameters, the trends of J as function 
of the DQD electrostatic potential and magnetic field, 
which we believe is helpful for the understanding and 
design of such devices. 

The results we present can only be treated semi- 
quantitatively for several reasons. First, the exact form 
of the potential is unknown and the problem is only ap- 
proximately two-dimensional. In an actual device, the 
dots will not be perfectly parabolic and the electrons are 
confined to a finite width in the direction perpendicular 
to the surface. Second, we use nine Gaussian basis el- 
ements per dot, and even though this results in energy 
levels converged to approximately 0.5% of their value, 
the exchange energy can have substantial fractional er- 
ror when it's absolute value is small (below 100 /ieV for 
our purposes). Even so, our results give a more accurate 
qualitative and semi-quantitative picture than previous 
variational approaches, and are sufficient to resolve the 
features of interest. In cases where additional accuracy 
is required, more powerful CI and iterative techniques 
(e.g. Poisson-Schrddingcr solvers) can be used. 9,10 Such 
techniques, however, are more computationally intensive 
and are best used to model realistic device structures 
rather than provide qualitative insight. 

While variations of the CI method with Hartree-Fock 9 
or molecular orbitals 11-14 have been used in the past, to 
our knowledge this is the first application of a CI method 
to a biased DQD capable of modeling both weakly and 
strongly coupled regimes. (The CI method in Ref. 14 is 
used to analyze an unbiased DQD, and is also restricted 
to the weakly coupled case.) 



FIG. 2: Spatial arrangements of Gaussian basis elements in a 
single quantum dot. Each solid black circle represents at least 
one basis element. There can be more than one element at a 
given point (circle) as specified by the table in Fig. 1. Note 
that these diagrams only describe the relative arrangement 
of the Gaussian centers; the "grid spacing" of the centers is 
optimized. 



A. Model 

A DQD qubit is considered within the effective mass 
approximation, and the electrostatic potential is approx- 
imated as the minimum of two parabolic dots, so that V 
has the form 



a wide range of physical systems, many of the solid-state 
approaches (e.g. Loss-DiVincenzo, 7 Kane, 8 and singlet- 
triplet 2 ), use a tunable exchange interaction between lo- 
calized electrons to perform qubit operations. The ex- 
change interaction causes a splitting between quantum 
states (with different spin) called the exchange energy, 
which we denote J. 

We focus here on a singlet-triplet qubit where two elec- 



V(x, y,z) = \ \m* x w\ min ((x — L) 2 + e, (x + L) 2 ) 

(5) 

+m* y LU 2 y y 2 + m* z u; 2 z 2 ] . 

This potential is parametrized in terms of e, L, and Q. It 
is readily seen from Eq. 5 that e is the energy difference 
between the left and right minima (the inter-dot bias), L 
is half the distance between the minima, and Eq = Huj is 
the confinement energy of the dots along the coordinate 



4 



axes. In the two-dimensional case, the potential is given 
by Eq. 5 with z set to zero. Cuts of the potential along 
the x-axis for different e arc shown in Fig. 3. We use 
GaAs material parameters, m* = 0.067 m e and k = 12.9 
unless otherwise noted. 
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FIG. 3: DQD potential along x-axis, V(x,0), for L — 30 nm 
and Eo = 3.0 meV. V is the minimum of parabolas centered 
at x — ±L with curvature proportional to Eq. 



In a real device, fabrication technique and gate volt- 
ages can effectively modify each of the parameters e, L, 
and Eq; the device design and voltage sequence will de- 
termine exactly how these parameters vary during device 
operation. 



B. Exchange Energy Results 

1. General Features 

Here we report, using GaAs material paramters, the 
general trends of J as function of the DQD electrostatic 
potential and magnetic field. We have in mind a system 
that varies e to perform qubit rotations, and so focus on 
the e dependence of J. Qualitatively similar behavior is 
seen in the Si/Si02 system, as will be shown in section 
IV B 3. 

First we consider the case of zero magnetic field. The 
shape of a typical J vs. e curve when B = is shown in 
Fig. 4 (for L = 30 nm and Eq = 3 meV) . The e-axis can 
be divided into three regions (marked in Fig. 4) corre- 
sponding to different charge character of the singlet and 
triplet states: (I) low e, where both singlet and triplet are 
in the (1,1) charge sector, (II) intermediate e, where the 
singlet takes on (0,2) character and the triplet remains 
a (1,1) state, and (III) large e, where both states are 
in the (0,2) charge sector. In (I), J is due to the differ- 
ence in Coulomb energy between spatially symmetric and 
antisymmetric compositions of single-dot wavefunctions, 
and scales exponentially with the inter-dot spacing (see 
Fig. 5). Note that at e — 0, we have dJ/de = exactly, 
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FIG. 4: Typical zero magnetic field exchange energy (solid 
line) as a function of dot bias. Dashed lines show the electron 
occupation of the left dot for the singlet and triplet state. The 
domain is divided into three intervals (I, II, III) based on the 
the charge configuration of the singlet and triplet states (see 
text). Electrostatic parameters are L = 30 nm, Eq = 3meV. 



due to the symmetry in e. In (II), J increases rapidly be- 
cause only the triplet is penalized for having an electron 
in the left dot. The curvature at the beginning and end of 
the linear portion is due to the inter-dot tunnel coupling 
(which increases with decreasing Eq and/or L). Larger 
tunnel coupling results in a more gradual transition be- 
tween (If) and the other two regions, as seen in Fig. 5. 
The exchange energy levels off in (III) since e does not 
affect the shape of the dot holding both of the electrons. 
The derivative dJ/de can be made arbitrarily close to zero 
by increasing epsilon, and in Fig. 4 decreases to O(10 -6 ) 
by e = 12 meV. The nearly constant value of J in (III) is 
essentially the exchange energy of a doubly-occupied sin- 
gle dot with confinement energy Eq. Figure 5 shows that 
J is insensitive to L in this region (since L does not affect 
the dot shape), and that increasing Eq raises the value 
of J. Note also that the boundary between regions (II) 
and (III) moves to larger e as E is increases, so that at 
large enough Eq the exchange energy actually decreases 
as Eq increases (see Fig. 5c). 

Now consider the case of finite magnetic field. In a 
single dot with two electrons, J oscillates as a function 
of B. 15,16 This is due to a competition between Coulomb 
and rotational energies which favors a state of increas- 
ingly higher angular momentum I. The ground state 
therefore alternates between an even-/ singlet and an odd- 
l triplet state, resulting in an oscillatory J. Increasing B 
also increasingly confines the electrons and contributes 
to a reduction of | J|. For a DQD with one electron in 
each dot, J is also a decaying oscillatory function of B 
which can be explained similarly. This behavior is shown 
in Fig. 6a for fixed Eq and L, and for e at representative 
values in intervals (I), (II), and (III) of Fig. 4. We see 
that for a fixed B, the exchange energy can have three 
qualitatively different types of dependence on e: (i) J 
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increases, (ii) J decreases then increases, and (iii) J de- 
creases. These three scenarios divide the B-axis into the 
intervals (i)-(iii) shown in Fig. 6a. Curves of J vs. e for 
B in each of these intervals are shown in Fig. 6b. Most 
noteworthy is case (ii) , in which J possesses a local mini- 
mum. This minimum exists because the B-field required 
to push J < is less in the (1,1) than (0,2) charge sector, 
but the derivative \dJ/dB\ is larger in the (0,2) sector. 10 



2. Convergence 

As the number of basis elements increases, the eigenen- 
ergies converge to the exact spectrum of the many-body 
Hamiltonian. Because of the basis optimization that is 
performed, a relatively small number of Gaussian basis 
elements (n g ps 20) result in energies converged on the 
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FIG. 6: Exchange energy as a function of magnetic field for 
L = 30 nm and Eo — 3meV. In (a), e = 0,3, 5meV corre- 
spond to regions I, II, and III of J vs. e curves respectively. 
The three different sequences of the curves in (a), delimited 
by the intervals (i), (ii), and (iii), lead to three qualitatively 
different types of J vs. e behavior. This is shown in (b), where 
5 = 1, 1.4, and 2T fall in (i), (ii), and (iii). The inset of (a) 
is a magnification around interval (ii). 



scale of several micro electronvolts. Figure 7 shows an ex- 
change vs. bias curve using n g = 18, 26, and 29 basis el- 
ements. The arrangements of these elements correspond 
to Fig. 2b on each dot, Fig. 2c on each dot, and Fig. 2c 
on each dot with 3 elements between the dots, respec- 
tively. As the insets of this figure reveal, the low-e region 
converges more rapidly than the high-e region. This is 
expected since in the low-e region the electrons are in sep- 
arate dots and the multi-electron wavefunction is well ap- 
proximated by linear combinations of the single-particle 
states in each dot, for which the Gaussian basis elements 
are excellent approximations. The high-e regime is that 
of a single doubly-occupied dot, which was considered 
in section III. We note that even though the energies 
have only converged to tens of /icV, the exchange energy 
has converged to approximately 5 /xeV due to error can- 
cellation. Since the energies have magnitudes of order 
lOmeV, the relative error is AE/E ps 0.05%. 



3. Differences in Si/SiCh 



FIG. 5: Zero-field (B — 0) behavior of the exchange energy 
when varying inter-dot separation (2L) and confinement en- 
ergy (Eo). In (a) J vs. e curves become sharper steps and J 
in the low-e region (I) decreases as L increases. This is illus- 
trated by J vs. L curves in (b), which decrease at low e and 
become flat at larger e, when the DQD is in the (0,2) charge 
sector. In (c) it can be seen that with increasing Eo, J vs. e 
curves rise at larger e and flatten out at larger J. Plots of J 
vs. Eo, shown in (d), increase until Eo the DQD transitions to 
the (1,1) charge sector. In (a) and (b), Eo is fixed at 3meV, 
and in (c) and (d), L — 30. 



Without any modification to our method, results for 
a Si/Si02 system in the single-valley approximation 
can be obtained using appropriate material parameters 
((m*,m*,m*) = (0.19, 0.19, 0.98) m e and re = (k Si + 
KSi0 2 )/2 = 8.0, which accounts for the image charge 
in the oxide 17 ). For the single- valley approximation to 
be valid, the splitting between states with different val- 
ley character must be large compared to the splitting 
between the ground and excited states of interest. We 
note that because the definition of Eo depends on the 
effective mass, the same quantum dot potential speci- 
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FIG. 7: Convergence of the CI method shown by comparing 
the exchange energy given by Gaussian bases with na ele- 
ments. Parameters Eo = 3meV, L = 30 nm, uj x = u) y , and 
B — are used, and represent the typical convergence behav- 
ior seen at other parameter values. Insets show zoomed views 
of the low- and high-e regions of the main plot. 
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FIG. 8: Comparison of the exchange energy obtained using 
material parameters m* and k appropriate for GaAs and Si02 
systems. DQD potentials with three different curvatures CI, 
C2, and C3, are used to generate the three curves for each 
material. These curvatures correspond to Eo = 1, 3, and 
5 meV in the GaAs system, but Eo will be different for the 
S1O2 system as described in the text. All DQD potentials use 
L = 30 nm and B = 0. 



fied by EQ aAs in a GaAs system will be given by Eq 1 
in a silicon system, where E^ = E§^ A8 \J mf aAs / mf* . 
Thus, the corresponding E values for E — 1, 3, and 
5meV in GaAs are E = 0.594, 1.782, and 2.969 meV in 
a Si/Si02 system. As expected, the exchange energy of 
an Si02 system shows qualitative behavior identical to 
GaAs systems. Quantitatively, the larger effective mass 
and smaller dielectric constant together result in lower 
values of the exchange energy given the same DQD po- 
tential (see Fig. 8). This comparison is somewhat arti- 
ficial, however, since different devices and experimental 
parameters (e.g. gate voltages) would be required to cre- 
ate identical potentials in the different materials. 



Figure 9 shows a comparison of the CI, HM and HL 
methods for the parameters of Fig. 4. Neither the HL 
nor HM method can predict the large-e flat since neither 
include triplet (0,2) states. While the HM method cap- 
tures the sharp rise in J (since it includes a (0,2) singlet), 
the onset of the rise is at much larger e. This is due to the 
tunnel barrier between the dots begin effectively larger 
because of the absence of basis elements which have sup- 
port between the dots. Said another way, in the HM case 
an electron can be centered in either the left or right dot 
but not between them, and therefore it requires a larger 
bias to push the electron out of the right dot. 
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4- Comparison with other approaches 

Lastly, we compare the CI method with Heitler Lon- 
don (HL) and Hund Mulliken (HM) techniques which 
previous studies 17-23 use to study DQD exchange en- 
ergy. The CI outlined here is more general than these 
methods from a variational standpoint: The full CI is a 
variational method and the space of trial wavefunctions 
includes the HL and HM spaces as long as the number 
of Gaussian functions uq > 2 (and is identical to HM 
for uq = 2). Effective Hubbard models are usually less 
precise than HM, since they approximate the Coulomb 
interaction as short-ranged. The small variational spaces 
of these methods restrict their applicability to qualitative 
results for weakly coupled quantum dots, 24,25 though the 
addition of variational parameters to the basis improves 
their accuracy. 15 ' 26,27 
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FIG. 9: Comparison of configuration interaction (CI), Hund- 
Mulliken (HM), and Heitler-London(HL) methods for E = 
3meV, L = 30 nm, and B = 0. Neither HL or HM methods 
can capture the flattening of the J vs. e curve at large e. 
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V. CONCLUSION 

Wc have described a configuration interaction method 
which uses Gaussian functions as basis elements and ap- 
plied it to two-electron single- and double-quantum-dot 
(DQD) systems, the latter functioning as a singlet-triplet 
qubit. The increased accuracy and extended regimes of 
validity of the CI relative to HL and HM methods, cou- 
pled with its ability to explore parameter space much 
more rapidly than mesh-based solvers, make the method 
uniquely suited to analyze the exchange energy in DQDs. 
We explain the general behavior of the exchange energy 
as a function of DQD parameters which correspond to 
inter-dot bias, inter-dot separation, dot size, and mag- 
netic field. In particular, we are able to analyze the 
transition from (1,1) to (0,2) charge sector, where more 
approximate methods fail, for both weakly and strongly 



coupled dots. This reveals new qualitative features of 
the exchange curve which could be useful for suppress- 
ing the effects of charge noise in physical qubits. 1 We 
focus on two-dimensional dots in a GaAs system, and by 
considering Si02 systems find that the qualitative fea- 
tures of the exchange energy are insensitive to typical 
changes in material parameters. The semi-quantitative 
results given in this analysis can be used to guide the 
tuning of experimental DQD devices seeking to imple- 
ment a singlet-triplet qubit . 

This work was supported by the Laboratory Directed 
Research and Development program at Sandia National 
Laboratories. Sandia is a multiprogram laboratory op- 
erated by Sandia Corporation, a Lockheed Martin Com- 
pany, for the United States Department of Encrgys Na- 
tional Nuclear Security Administration under Contract 
DE-AC04-94AL85000. 



Appendix A: Matrix Elements between s-type Gaussian functions 

1. IP Gaussian matrix elements 

We derive here the matrix elements for single-particle operators which are piecewise polynomial in the components 
of particle position r and momentum p. Since the basis functions are Gaussian, this reduces to the case of operators 
piecewise polynomial in the components of f (p = —iKSJ brings down powers of r from a Gaussian's exponent). Let 
us define two Gaussian basis elements, \gi) and \gj), by their real-space representations, 

g i {f)=Ne-^- ?A ^ ( - p - fA) , gj{r) =N'e-^- fB ^ i> - i>B) (Al) 

Vectors ta and fg are the center positions of \gi) and \gj), and a and f3 are d x d diagonal matrices specifying their 
exponential factors. We will make repeated use of the identity 

9i(r)9j(r) = NN' 'Jjr e -(*-Ai*W-Au») ( A2 ) 

where 

K = e -l?A-?B)c(f A -i'B) with C = a/3/(a + /3) (A3) 

Rab = {ar A +Pr B )/{a + P) (A4) 

ft = a + fi (A5) 

which allows us to transform the product of two Gaussians into a single Gaussian. Throughout this appendix, division 
by a matrix means multiplication by its inverse. 

Next, let us define a "piecewise-polynomial" operator O as one which can be written in the form 

O = Y^Oitfxfr.wW (A6) 

i 

= £ £*IK M) (at) 

i \ t k I 

Oi is polynomial in the components of r, and is expanded in polynomial terms in the second line. The charac- 
teristic function X\st\ equals 1 within the d-dimensional interval [a,b] (e.g. a cube for d = 3) and everywhere 

else. The components of a and b are real or ±oo. Thus, (gi\0\gj) is a linear combination of terms with the form 
(5i|X[s 3 (?) Ilfe=i r k k \9j) where n, a, and b are d-dimensional vectors. Wc compute this general element immediately 
and then use the result to find the matrix elements of the kinetic and potential energy operators used for DQD 
Hamiltonians. 



8 



We begin the computation as follows 

d 



tap^l f[ r k k \9i) = NN ' K I d d ue~^ a f[ (u k + R ABk ) n » 
k=i Js ' fc=i 



NN'K U f" du k e-^< (u k + R ABk ) nk 



d n 



"fc / \ rb, 
NN'KY[Y^ ( ) {RABkT k ~ l I 



du k u k e 



l „-fJ-kk'U. k 



k=l l=Q 
d n 



k=l 1=0 



NN'K [] E ( ?' ) (RAB k ) nk ~ l F(l; a' k , b' k ,» kk ) 



(A8) 
(A9) 
(A10) 
(All) 



In the first line we have used (A2) and made the substitution u = r — Rab- The limits of integration are from 
a' = a — Rab to b' = b — Rab- hi the second line the integral written as product of one-dimensional integrals, and 
in the third line the binomial expansion of (u,; + RABi) ni is inserted. The fourth line (Eq. All) defines the integral 



F(l;a,b,u) = / duu l e-^ 

J a 



(A12) 



which we now compute. 

Using the shorthand notation F(l) for F(l; a, b, /i),and integrating by parts {U = dV = ue^ u du, dll 

{I - l)u l ~ 2 du, V = -l/(2/^)e-^ 2 ) 



F(l) 



2/i 



a + 2/X J a 



i^e'^du 



,1-1 



/ - 1 — w 
G{1) + — — F(l - 2) where G(l) = 



i b 



2/i 



Since 



F(0) 
F(l) 



/ e ^ du = /— [crf(^6) - cr^^a)] and 

J a * V L 1 



ue~^ u du 



1 r 

271 



can be computed directly, we can write F(l) non-recursively by 



i/2 



Z - 2p + 1 



P =i 



2/i 



(7 even) 



where 



Z/2 — 1 Z even , F* — I ' even 



(Z - l)/2 — 1 I odd 



F(l) Zodd 



(A13) 
(A14) 

(A15) 
(A16) 

(A17) 

(A18) 



In practice, F(l) is calculated by a routine which iteratively builds the solution, such as the following (pseudo 
C++): 

real F(k) { 

i = ( k is even ) ? : 1 
x=F(i) 



while ( i < k ) { 
i++ 
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x *= i/(2*mu) 
i++ 

x += G(i) 

} 

return x 



With the above analytic expression of F(l; a, b, fi), Eq. All can be used to compute (gi\X[z g] (^0 r k k ano - ^ nus 

the matrix elements of any "piecewise-polynomial" operator O. 

Note that if the matrix element is a single polynomial, so that a k = — oo and b k = +00 for k = 1 . . . d, then F(l) is 
zero for odd I (integrand is odd) and for even I is given by the elementary integral 



x 2n e -ax- 



dx = 2vtt — r 



1 



2^fd 



2)1+1 



(A19) 



Thus, in this simplified case when there are no characteristic functions. 



(^1 n 



r T \9j 



= NN'K \ 



k=l 



k=l 



/— O.cvcn 



{RABkf 



II 



1 



(1/2)1 \2^jm 



1+1^ 



(A20) 



Piecewise-polynomial potentials, such as parabolic or quartic quantum dot potentials, are naturally expressed in 
the form given by Eq. A6. Hamiltonian terms involving the momentum operator, however, require a few preliminary 
steps. For example, matrix elements of the kinetic energy operator, (gi \ — V 2 |gj), can written as matrix elements of 
the operator (r — r A )(Aa(3)(r — r B ), which is polynomial in r by integrating by parts and taking derivatives of the 
Gaussian basis functions. 



V 2 \ 9j ) = J d d rV ai (r)-V 9l {r) 



(A21) 

= / d d r gi {r)g {r){r-r A ){AaP){r-r B ) (A22) 
= { gi \{r ~ ?A)^a(i)(r - r B )\ 9] ) (A23) 

Next define diagonal matrix P = 5/3, and use Eq. A20 to arrive at a formula in terms of the basis element parameters 

(A24) 

(A25) 



(9i\ ~ V%j) = 4(g t \rPr- fPf B - f A Pf + f A Pr B \ gj 

d 



Pkk R ABk 



_?L + _v!L TT — 



[r Ak P kk r Bk - Pkkr B k - r Ak P kk ) TT J — 

, V Mi/ 



ic Pkk 



2 f i 



h II \/J + ( R ABk - (r-Ak + r Bk ) + r Ak r Bk ) J] J^L 



kk l=£k 



(A26) 
(A27) 
(A28) 



We define K, = ANN'K in Eq. A27, and in the last line we have defined A A = R AB — r A and A B = R AB — r B . 
The kinetic matrix element may also be computed directly as follows: 



(.9)1 -V 2 | 5j ) 
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d d fV gi (r) • V.9j(f) 

(f - itts) + (R AB - Fa)] P \(F- R AB ) + (Rab - P B ) 



K J d d re- (f - RAB)Jl 
K I d d ue~ ajlS 
K 



uPu + uP^ — -(Fa -Fb) + (Rab - Fa)P(Rab - Fb] 

a + /3 



K, 



K. 



K. 



J2 J d d uPuU 2 ie~^ + A A PA B J d d ue 

.1—1 

(J du^e-^A (n/ " "') +A A PA B f[ J d Ui e-^< 
i= i \ 2 Mii / j^i v ^ i= i V M« 



(A29) 
(A30) 
(A31) 

(A32) 
(A33) 
(A34) 

(A35) 
(A36) 



We integrate by parts to get (A29), use (A2) to transition to (A31), and substitute u = r — Rab in (A32). The term 
of (A32) linear in u vanishes since the integrand is odd. We then expand the matrix notation to arrive at (A35). The 
final line is obtained using the elementary integrals 



e ax dx 

CO 

x 2 e ax2 dx 



7T 

a 
n 



2a 3 / 2 



(A37) 
(A38) 



The overlap matrix element can be calculated using Eq. All or A20, though it is straightforward to compute 
directly by using Eq. A2 and the translation invariance of the integral. Following the latter approach, we obtain 



(9i\9j) = NN ' / d d rKe- r ^ r = NN'Kt 



det [i 



(A39) 



2. Coulomb matrix elements 

To compute the matrix elements of the m-body Hamiltonian (1) in the basis B^ B , matrix elements for the two- 
particle Coulomb term — must be computed. We introduce two more Gaussian basis elements, \gi>) and \gy) written 
similar to those in (Al): 

gi ,(r) = N"e- {p - pc)W - pc) , gj,(F) = }?» e -<?-*i>W-*D) (A40) 
Anticipating the combination of <7i'(f) and gj'(r) using (A2), define the following analogous to Eqs. (A4)-(A5): 



K' = e ~( p o-ro)C(rc-P D ) with c = 78 / ft + 8) 
Rcd = {ir c + 5r D ) I (7 + 5) 
v = 7 + S 



(A41) 
(A42) 
(A43) 



We now turn to the matrix element of interest. We begin by writing it as a real-space integral then apply Eq. A2 to 
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each pair of Gaussian basis elements: 



(9i9i'\ — \9j9f) = — d rid f 2ffi (fi)5i/(f2) — gj{ri)gj>{r 2 ) 



1 



]Q_ I d d rid d r2e~^ 1 ~ RAB * > ^ Fl ~ RAB * > e ~(r2-RAB)fj.(r 2 -RAB) 



(A44) 
(A45) 



where JC = NN' N" N'" KK' . Next, we write each of the two exponentials and \jr\% in terms of their Fourier 

transforms, defined by FT [/(f)] = (2ir)~ d f d d kf(k)e lk ' r . At this point we fix d = 3, so that the Fourier transform 
of 1/r is well defined. The 2D solution will be obtained later, by taking a limit of the 3D result. Thus, the Fourier 
transforms 



FT 



det i 



e 4 



and 



FT 



47T 

— for a = 3 

k l 



are inserted into in Eq. A45 and result in (for d = 3) 

e 2 JCe 2 

{gi9i'\—\gjg 3 ') = , x 3rf 



d d r 1 d d r 2 d d k 1 d d k 2 d d k 3 



71 -ifci^r'fci 



det /j, 



-,e 4 



47T / 7T 



x ,/ _ e ~\k 3 v 1 k3 e ik 1 -(r 1 -R AB ) e ik2-(f 1 -r 2 ) e ik3-(r2-R.CD) 



k 2 V detf 



(A46) 
(A47) 

(A48) 
(A49) 



Integrating over fi and f% yields delta functions (2-7r) d (5(fci + k 2 ) and (2ir) d 5(k 2 — k^). Then integrating over k 2 and 
&3 effectively sets —k\ = k 2 = k% in the integrand, and we define k = —k± to clean up the notation. After these 
integrations, Eq. A49 becomes 



(9i9i'\ — \9j9j') - 

kt (2,7rj k det fit 



d dk^- e -fck e ik.& 



(A50) 



where a = ^tif (division is multiplication by matrix inverse) and A = Rab — Rcd ■ Define a = Tr(a)/d, write 



e = exp(— ak ) exp —k(a — al)k (1 is the identity matrix in d dimensions), and use the identity 

2a / dSS~ 3 e-° k / s 



k 2 

(this follows from e~ a = a e~ ax dx with x — l/S 2 and a = ak 2 ) to transform Eq. A50 into 



(9i9i'\ — \9j9j' 
nr 



JCe 2 4tt 7r fl 



2a I 1 S / d d ke- k ^+^ s ^ k e lk - 1 



« W^det^ 7o S 3 
JCe 2 Air n d ^ f 1 d^ n^e^f^+^^ 
« (27r) d VdeTp Jo S 3 V det ( ? - (J + cr / 5 ' 2 ) 

JCe 2 



hi 



iT 5/2 a (det/iv) 2 I (a, A) 



(A51) 

(A52) 
(A53) 

(A54) 
(A55) 



where <r = crl when used in a matrix context, and in the last line we explicitly put d = 3. In Eq. A54 we have defined 
the integral 



7(5, A) 



1 dS_ 

o 

1 dS_ 

o 



[det(5 - a + a/S 2 )p g-^C^-"^/^)- 1 * 



n(^- CT (l-V5 2 )) 



-iEI^Af/Ca.i-o-a-l/S 2 )) 



(A56) 
(A57) 
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The second line follows since a is a 3 x 3 diagonal matrix. We cannot express the integral in closed form, and so must 
compute I(a, A) numerically. In the two-dimensional case, we take the limit /Z33 = 1/33 — > 00, which means that 033 — > 

0. This is the limit where all the Gaussian basis elements have the same width 77 in the z-dircction, and 77 approaches 
zero (making the elements two-dimensional). K. contains the factor (2a3 3 /n) 1 / 4 (2l3 33 /'K) 1 / 4 (2j 33 /T:) 1 ^ 4 (2S 33 /7:) 1 ^ 4: = 
277/77, which cancels the factor 1/^/^33^33 = l/V (277) (2?;) = 1/(277) from (det jlv^ 1 / 2 , leaving a factor of 7T in the 
denominator. This reduces the exponent of 7r in Eq. A55 from 5/2 to 3/2 as seen below in Eq. A59. 

The integral I(a, A) converges. The only possible trouble occurs when (an — a(l — 1/S 2 )) = 0, or equivalently 
5=1/ s/T— an/a, for some i = 1, 2, 3. Since an > implies that 1/ y/l — an /a is either imaginary or greater than 

1, the only divergence of the integrand for S € [0, 1] occurs at S = 1 when an = 0. Although this happens for i = 3 
in the 2D case (0-33 = 0), the divergence is ~ l/\/l — S, which is intcgrablc. Thus, the integral is well defined and 
convergent over the entire range of physical parameters. 

In summary, the Coulomb term matrix elements for two and three dimensions can be carried out analytically up 
to the numerical evaluation of a convergent one-dimensional integral. They are given by: 



\9i9i' —\9j9j> 



(9i9i'\ — \9j9j' 

fxl 



— Tr 5/2 a {det fiVy * I(a, A) (3D) 

K 



K-2De 2 3/2 



a ( det jjiv 



/(a, A) (2D) 



(A58) 
(A59) 



The subscripts "2D" are a reminder that the normalization factors and determinant contain only x and y factors. 



Implementation Note: 

We note additionally that when the integral for two dimensions is computed numerically, the integral around 5=1 
is approximated by a closed form. At S = 1 — e, when e -C 1, the integrand in Eq. A56, which we denote W(S), is 
approximated via the expansion 



a-a\\- 



2ae 



(A60) 



resulting in an approximation for the integral between S = 1 — £q an d 1. 



dSW(S) 



dS det {a + a{2e)y 1/2 (1 - e )- 3 e -i( A ^ 1+A >v ') 

— e 

£ ° de (a x a y a(2e)y 1/2 e ~i( A >* 1+A >^ 
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e~4 



HK^+KO r° de 



y / 2a x a y a 



^Ja x a y a 



(A61) 
(A62) 
(A63) 
(A64) 



3. Isotropic Gaussian functions - the completely analytic case 



When all the Gaussian basis elements arc isotropic, that is, when 5 is a proportional to the identity, then Coulomb 
matrix elements can be computed analytically in both two and three dimensions. In 3D, we follow the above derivation 
up to Eq. A50 and let ju = /jl, v = ul, and a = al define the scalars 71, v, and a corresponding to the similarly 
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named matrix. Then write Eq. A50 in spherical coordinates with the ^-axis along A to obtain 



(9i9i'\ — \9j9j') = fC' I k 2 dkd(cost 



4tt 



^ — ak 2 gifcA cos 6 



= AmKJ \ dk 



-ah 



poo p -ak~ 

8ttJC' / dk 



ikA 



d(cos0) e 



ikA cos 6 



2,T 



= IGjtK' / dke~ ak2 



ikA 



f A f 00 i 
16tt/C' dy dke- ak cos ky 
Jo Jo 

f A f°° 2 

8nJC' dy dke- ak cos ky 

Jo J -oo 



(A65) 
(A66) 

(A67) 
(A68) 
(A69) 
(A70) 



where K! = ^(^)" 3/2 . In line A69 we introduce a dummy variable y in order to make the fc-integral tractable (see 
below). The final line uses the fact that the integrand is even to extend the range of integration. By writing cosfcy 
in exponential form and completing the squares, we can integrate over k to get 



{9i9v\^ r \9i9i') = 47r ^' / dv I dke-° k2 (e ik v + e- ik «) 



4tt/C' / dy 



dk 



/•A poo 2 

4tt/C' dy dfc2e- CTfe Vfe 

Jo J -oo 



87T/C' 



8ttK\ - 



dye' 



crf(A) 



,2tt 2 
AC' erf(A) 

(7 



(A71) 
(A72) 
(A73) 
(A74) 
(A75) 
(A76) 



We obtain line A73 by shifting k in each of the terms (in different directions), which does not alter the integral and 
introduce the error function erf(x) = ^= J Q X dx' e~ x on line A75. Expanding Kf gives us a final analytic formula for 
Coulomb matrix elements in three dimensions {d — 3), 



(9i9i'\^\9j9f) = 4^ (^3/2 ei 



(A77) 



In 2D, we proceed to Eq. A45 as above, but then insert 2D Fourier transforms instead of the 3D ones yielding 
Eq. A49. The Fourier transform of the Gaussian basis elements has the same form as the 3D case, but now 



FT 



= — for d = 2 
k 



(A78) 



Inserting these into to Eq. A45 (with d = 2) gives 

^4f~ I d^d^hd^hi Ke-^^ 
lufn J V /i 2 



(9i9i'\ — \9j9j' 
nr 



JCe 2 

(27 

2tt In 2 



x . ; g -k 3 /(iu) ^^-(tx-Rab) e ik 2 -(ri-r 2 ) e ik 3 -(r 2 -RcD) 

kl V V 2 



(A79) 
(A80) 
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As in the 3D case, integrating over T\ and f*2 produces delta functions which are removed by then integrating over ki 
and &3. This effectively sets k = —k\ = k^ = k% in the integrand, 



(9i9i'\ — \9j9j' 



ICe 2 1 
4k pM 



k 



(A81) 



where a = and A = Rab — Rcd- Again we use the scalars u, v, and a which correspond to matrices in the more 
general (non-isotropic) case. Changing to polar coordinates with the x-axis along A, 

„2 



(9i9^'\ — \9j9 3 ' 



ICe 2 


n 


2k 


//;/ 


ICe 2 


7T 


2k 


\XV 


ICe 2 


7T 


2k 


//;/ 


ICe 2 


7T 


2k 


//;/ 


ICe 2 


TT 2 


2k 


[IV 



(W 

de 
de 



dk i 



dk 



-<j( k- 



e- A l^e 



2tt 



dO'e- 



~A 2 /(8cr)j o 



-A 2 
~8cT 



(A82) 
(A83) 
(A84) 
(A85) 
(A86) 



On line A83 we have completed the square, and to obtain line A84 the trigonometric identity cos 2 6 = (1 + cos 29) /2 
is used. The quantity in square brackets on line A84 is equal to \tkJc! (the variable of integration can be shifted). 
The parenthesized quantity on line A85 is equal to J* d9' exp [— (A 2 /8ct) cos by symmetry, which is equal to 
ir Iq(— A 2 /(8a)) where Iq is the first modified Bessel function. Combining terms, we arrive at a final analytic formula 
for Coulomb matrix elements in 2D, 



(9i 



g v \-\9 3 9y) = ^—e^l^ I (-£-) (2D, isotropic) 

kt 2k^Jo [IV \ 8(7 J 



(A87) 



Appendix B: Optimization of Gaussian function parameters 



The computational speed of the configuration interaction method described in these notes is dependent on its 
ability to use relatively small basis sizes (i.e numbers of Gaussian elements) to obtain semi-quantitative results. To 
achieve a more rapid convergence with respect to the basis size, the placement and exponential factor of the Gaussian 
basis elements is optimized within a subspace of all possible sets. This optimization process is of great practical 
importance, since results must converge before reaching the maximum uq allowed by our computational resources 
(currently 50 — 100). In this section, we explain in detail the method used to choose an "optimal" set of Gaussian 
basis functions. 

The number of basis elements tiq is always given as a fixed parameter to the optimization procedure. An initial 
basis of nc elements is generated by specifying either (I) the location and size of each dot, along with and the number 
of elements to place in it, or (II) the location (in a 2D plane) and size of each Gaussian function. Method (I) requires 
that the system be comprised of one or more quantum dots, where method (II) can be used for any system. 

Positions of the elements in a Gaussian basis are generated from an underlying two-dimensional mesh of points with 
two characteristic length scales a x and a y (usually the spacing between elements associated with the same quantum 
dot along the x- and y-axis, respectively). Method (I) creates this mesh based on location and size of the dots, and 
always places a mesh point at the center of each dot. In method (II), the locations of elements are given in terms of 
a x and a y , which are given separately. 

The elements in a basis are partitioned into nss subsets such that the exponential factor, is the same for elements 
in the same subset. Initial values for these factors, denoted on for i = 1 . . .nss, are either chosen based on the dot 
size (method I) or specified directly (method II) . If the Gaussian elements are required to be isotropic, each a must 
be a multiple of the identity. This restriction results in less freedom for basis optimization but increased computation 
speed. In the case of method I, nss = 2, and the elements at the center of each dot are allowed to have a different 
exponential factor than the rest of the elements. When there are multiple basis elements at a given mesh point the 
exponential factors of additional elements are found by multiplying the previous element's coefficient by a constant 
factor A. Method I fixes A = 1.5 and method II takes a value for A as input. 
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Optimization of the basis is performed by minimizing an energy E with respect to a x , a y , a,, and A simultaneously. 
The energy minimized is the either (1) the lowest single-electron energy, (2) the lowest many-body energy, or (3) the 
lowest many-body energy with a given symmetry (e.g. S z = and total spin S = 1). In the results for the exchange 
energy presented in this work we perform minization via case (3) twice: once to minimize the lowest singlet energy 
and once to minimize the lowest unpolarized triplet energy. 
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